Нам понадобятся следующие пакеты.

1. Описание датасета:

1.1. Всего 72 мыши, у каждой по 15 измерений.

1.2. Можно выделить 8 классов, которые различаются по генотипу (контроль, трисомия), лекарство (memantine (m) or saline (s)), поведение (context-shock (CS) or shock-context (SC)). Итого 8 классов: c-CS-s, c-CS-m, c-SC-s, c-SC-m, t-CS-s, t-CS-m, t-SC-s, t-SC-m.

1.3. Мыши по группам не сбалансированы: в группах от 7 до 10 мышей.

c-CS-s: 9 mice c-CS-m: 10 mice c-SC-s: 9 mice c-SC-m: 10 mice

t-CS-s: 7 mice t-CS-m: 9 mice t-SC-s: 9 mice t-SC-m: 9 mice

1.4. Количество полных наблюдений:

d <- d[,-1]

aggr(d, sortVars=F, combined=T, bars=F, numbers=T, prop=F, sortCombs=T) #see NAs

Мы видим 3 наблюдения с NA в больше половине случаев, их я удалю. Всего 552 наблюдения с полной информацией.

d <- d[- (which(rowSums(is.na(d))>10)),] ##Exclude obs with many NAs
print("Количество оставшихся NA")

## [1] "Количество оставшихся NA"

sort(colSums(is.na(d)),decreasing=T)[1:10]

##     BCL2_N   H3MeK4_N      BAD_N     EGR1_N  H3AcK18_N    pCFOS_N      ELK_N 
##        285        270        213        210        180         75         15 
## Bcatenin_N      MEK_N   DYRK1A_N 
##         15          4          0

2. Различия в уровне продукции BDNF_N в зависимости от класса.

Посмотрим на pointrange.

d_d <- d
d_d$class <- as.factor(d_d$class)
ggplot(d_d, aes(x = d_d$class, y = d_d$BDNF_N, colour = class)) +
  labs(y = 'BDNF_N', x = 'Class') +
    theme_light() +
  stat_summary(geom = 'pointrange', fun.data = mean_cl_normal)

По графику видно, что группы из середины отличаются. Сделаем дисперсионный анализ

# функция Anova из пакета car
mod_1 <- lm(BDNF_N ~ class, data = d_d)
m_anova <- Anova(mod_1)

# встроенная функция
m1 <- aov(d$BDNF_N ~ d$class, data = d_d)
summary(m1)

##               Df Sum Sq Mean Sq F value Pr(>F)    
## d$class        7 0.2878 0.04112   18.82 <2e-16 ***
## Residuals   1069 2.3362 0.00219                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# просто сравнила функции - показывают одно и тоже :)

По данным дисперсионного анализа, мы видим, что уровень BDNF_N достоверно зависит от класса животного.

# Данные для графиков остатков
mod_diag <- fortify(mod_1)
# График расстояний Кука
ggplot(mod_diag, aes(x = 1:nrow(mod_diag), y = .cooksd)) +
geom_bar(stat = 'identity')

Остатки распределнены нормально.

ggplot(mod_diag, aes(x = mod_diag$class, y = .stdresid)) + geom_boxplot()

Остатки по классам не различаются.

qqPlot(mod_1, id = FALSE)

Распределение очень близко к нормальному.

Сделаем пост-хок тест, чтобы узнать, какие группы различаются. Я взяла тест Тюкей, поскольку он относительно не строгий и не слабый.

posthoc <- glht(mod_1, linfct = mcp(class = "Tukey"))
summary(posthoc)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = BDNF_N ~ class, data = d_d)
## 
## Linear Hypotheses:
##                        Estimate Std. Error t value Pr(>|t|)    
## c-CS-s - c-CS-m == 0  0.0030979  0.0055459   0.559  0.99930    
## c-SC-m - c-CS-m == 0 -0.0482717  0.0053980  -8.942  < 0.001 ***
## c-SC-s - c-CS-m == 0 -0.0258249  0.0055459  -4.657  < 0.001 ***
## t-CS-m - c-CS-m == 0 -0.0264852  0.0055459  -4.776  < 0.001 ***
## t-CS-s - c-CS-m == 0 -0.0337570  0.0059483  -5.675  < 0.001 ***
## t-SC-m - c-CS-m == 0 -0.0181541  0.0055459  -3.273  0.02392 *  
## t-SC-s - c-CS-m == 0 -0.0136310  0.0055790  -2.443  0.22105    
## c-SC-m - c-CS-s == 0 -0.0513696  0.0055459  -9.263  < 0.001 ***
## c-SC-s - c-CS-s == 0 -0.0289228  0.0056900  -5.083  < 0.001 ***
## t-CS-m - c-CS-s == 0 -0.0295831  0.0056900  -5.199  < 0.001 ***
## t-CS-s - c-CS-s == 0 -0.0368549  0.0060829  -6.059  < 0.001 ***
## t-SC-m - c-CS-s == 0 -0.0212520  0.0056900  -3.735  0.00474 ** 
## t-SC-s - c-CS-s == 0 -0.0167289  0.0057223  -2.923  0.06894 .  
## c-SC-s - c-SC-m == 0  0.0224468  0.0055459   4.047  0.00155 ** 
## t-CS-m - c-SC-m == 0  0.0217865  0.0055459   3.928  0.00226 ** 
## t-CS-s - c-SC-m == 0  0.0145147  0.0059483   2.440  0.22285    
## t-SC-m - c-SC-m == 0  0.0301176  0.0055459   5.431  < 0.001 ***
## t-SC-s - c-SC-m == 0  0.0346406  0.0055790   6.209  < 0.001 ***
## t-CS-m - c-SC-s == 0 -0.0006603  0.0056900  -0.116  1.00000    
## t-CS-s - c-SC-s == 0 -0.0079321  0.0060829  -1.304  0.89717    
## t-SC-m - c-SC-s == 0  0.0076708  0.0056900   1.348  0.87968    
## t-SC-s - c-SC-s == 0  0.0121939  0.0057223   2.131  0.39494    
## t-CS-s - t-CS-m == 0 -0.0072718  0.0060829  -1.195  0.93310    
## t-SC-m - t-CS-m == 0  0.0083311  0.0056900   1.464  0.82595    
## t-SC-s - t-CS-m == 0  0.0128542  0.0057223   2.246  0.32461    
## t-SC-m - t-CS-s == 0  0.0156029  0.0060829   2.565  0.16969    
## t-SC-s - t-CS-s == 0  0.0201260  0.0061130   3.292  0.02262 *  
## t-SC-s - t-SC-m == 0  0.0045231  0.0057223   0.790  0.99360    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Построим график, чтобы визуализировать результат (одинаковые буквы означают неразличающиеся группы)

MyData <- data.frame(class = factor(levels(d_d$class), levels = levels(d_d$class)))
MyData <- data.frame(MyData, predict(mod_1, newdata = MyData, interval = 'confidence'))


ggplot(data = MyData, aes(x = class, y = fit)) +
  geom_bar(stat = 'identity', aes(fill = class), width = 0.5) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.1)+
labs(x = 'Class', y = 'Экспрессия белка BDNF_N') +
geom_text(aes(y = 1.6, label = c('A', 'A', 'B', 'C', 'D', 'E', 'A', 'A')), colour = 'black', size = 5, position = position_stack(vjust = 0.5))+
  theme(legend.position = 'none')

3. Попробовать построить линейную модель, способную предсказать уровень продукции белка ERBB4_N на основании данных о других белках в эксперименте

Поскольку у нас около 75 белков, до построения модели, попробуем сократить количество анализируемых переменных. ### 1. Анализ на мультиколлинеарность. Попарно сравним коэффициент корреляции белков.

m <- length(colnames(d))
# Функция для расчета p-value для коэффициента корреляции (взято: http://www.sthda.com/english/wiki/visualize-correlation-matrix-using-correlogram)
# mat : is a matrix of data
# ... : further arguments to pass to the native R cor.test function
cor.mtest <- function(mat, ...) {
    mat <- as.matrix(mat)
    n <- ncol(mat)
    p.mat<- matrix(NA, n, n)
    diag(p.mat) <- 0
    for (i in 1:(n - 1)) {
        for (j in (i + 1):n) {
            tmp <- cor.test(mat[, i], mat[, j], ...)
            p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
        }
    }
  colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
  p.mat
}
# matrix of the p-value of the correlation
d <- as.data.frame(d)
p.mat <- cor.mtest(d[,-c((m-3):m)])
head(p.mat[, 1:5])

##              DYRK1A_N      ITSN1_N        BDNF_N         NR1_N        NR2A_N
## DYRK1A_N 0.000000e+00 0.000000e+00  3.436576e-34  5.382704e-23  6.337133e-28
## ITSN1_N  0.000000e+00 0.000000e+00  1.798784e-57  7.246317e-48  9.222316e-49
## BDNF_N   3.436576e-34 1.798784e-57  0.000000e+00 7.780844e-247 3.484596e-195
## NR1_N    5.382704e-23 7.246317e-48 7.780844e-247  0.000000e+00  0.000000e+00
## NR2A_N   6.337133e-28 9.222316e-49 3.484596e-195  0.000000e+00  0.000000e+00
## pAKT_N   2.215088e-09 1.111704e-06  1.177053e-26  2.320727e-12  2.918397e-04

cor <- cor(d[,-c((m-3):m)],use = "pairwise.complete.obs")


corrplot(cor, type = 'upper', p.mat = p.mat,
         sig.level = 0.01, insig = 'blank', 
         tl.col="black", tl.srt=45,
         tl.cex = 0.4)

Данных очень много, поэтому крайне высока вероятность наличия мультиколлинеарности. Одним из признаков мультиколлениарности является высокая парная зависимость между малозначимыми переменными, поэтому найдем пары белков с наибольшей коллинеарностью и удалим одного из пары.

## [1] "Number: Correlations between proteins > 0.8"

## [1] 87

## [1] "Number: Correlations between proteins > 0.5"

## [1] 435

87 белков из датасета имеют корреляцию больше 0,8. Удалим каждого из пары таких белков.

Var_na <- which(colSums(is.na(d))>0)
d <- as.data.frame(d)
cor_na <- cor(d[,-c(Var_na,(m-3):m)],d[,Var_na],use = "pairwise.complete.obs")
d <- d[, which(! colnames(d)  %in% names(which(colSums(cor_na > 0.8) > 0)))]

m <- length(colnames(d))
Vars_Mcol <- findCorrelation(cor(d[,-((m-3):m)],
                                 use = "pairwise.complete.obs"),0.8)
d <- d[,- Vars_Mcol]

Осталось 54 белков.

m <- length(colnames(d))
d <- as.data.frame(d)
aov <- numeric()
for( i in (1: (m-4))){
aov[i] <- summary(aov(d[,i]~d$ERBB4_N))[[1]]$"Pr(>F)"[1]
}
print("number of protein variables which have a significant relationship to ERBB4_N")

## [1] "number of protein variables which have a significant relationship to ERBB4_N"

sum( aov < 0.05/length(aov))

## [1] 37

Vars <- which(aov < 0.05/length(aov))
d <- d[,c(Vars,(m-3):m)]

Итого осталось 37 белков, которые значимо коррелируют с ERBB4_N. Для построения обычной линейной модели предикторов слишком много, сделаем РСА.

which(colnames(d)=='ERBB4_N')

## [1] 23

m <- length(colnames(d))
# уберем колонки c не белками
prot_only <- as.matrix(d[,c(-1,-23, -(m-3):-m)])
# стандартизируем переменные
# Function to standardize numeric variables
z_numeric_vars<-function(data_frame) {
  for (n in names(data_frame)) { 
    if (class(data_frame[[n]]) == "numeric" | class(data_frame[[n]]) == "integer"){
      var = paste(n,"_z", sep="")
      data_frame[[var]] <- scale(data_frame[[n]], center = TRUE, scale = TRUE)
      data_frame[[n]] = NULL
    }
  }
  data_frame
}
prot_st <- z_numeric_vars(prot_only)
# заменим NA на среднее
prot <- na.aggregate(prot_st)

Построим ординацию

Построим график собственных чисел

screeplot(ord, bstick = TRUE, type = 'lines')

## [1] 52.5942

Оставим первые 3 компоненты, они объясняют 53% изменчивости.

Построим факторные нагрузки

biplot(ord, scaling = 'species', correlation = TRUE,
main = 'PCA - species scaling', display = 'species')

Построим график для первых трех компонент.

df_scores <- data.frame(d, scores(ord, display = 'sites', choices = c(1, 2, 3), scaling = 'sites'))

plot_ly(df_scores, x = ~PC1, y = ~PC2, z = ~PC3, color = df_scores$class, size = 0.5)

prot_1 <- as.data.frame(prot)
ERBB4_N <- scale(d$ERBB4_N)
new_d <- as.data.frame(cbind(ERBB4_N, scores(ord, display = 'sites', choices = c(1, 2, 3), scaling = 'sites'))) 

lmod <- lm(V1 ~ ., data = new_d)
summary(lmod)

## 
## Call:
## lm(formula = V1 ~ ., data = new_d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.54187 -0.36359 -0.02816  0.34886  2.60972 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.829e-16  1.941e-02   0.000        1    
## PC1          3.365e+00  9.421e-02  35.716  < 2e-16 ***
## PC2         -7.180e-01  1.107e-01  -6.486 1.35e-10 ***
## PC3          2.133e+00  1.323e-01  16.125  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6371 on 1073 degrees of freedom
## Multiple R-squared:  0.5952, Adjusted R-squared:  0.5941 
## F-statistic: 525.9 on 3 and 1073 DF,  p-value: < 2.2e-16

Модель объясняет 59% изменчивости.

Проверка на мультиколлинеарность.

sqrt(vif(lmod))>2

##   PC1   PC2   PC3 
## FALSE FALSE FALSE

Мультиколлинеарности нет.

Однако, vif выдаёт такой странный результат. (если можно, поясни, пожалуйста, в каком месте я ошиблась или что с этим делать)

vif(lmod)

## PC1 PC2 PC3 
##   1   1   1

Проверим остатки на нормальность

qqPlot(lmod, labels = row.names(states), simulate = TRUE, main = 'График Q-Q')

## sit142 sit250 
##    142    250

Есть какие-то выбросы, но в целом похоже на нормальное распределение.

Проверим остатки на линейность.

crPlots((lmod))

Если нелинейные связи есть, то они не очень большие.

Построим визуализацию.

new_data <- data.frame(ERBB4_N = seq(min(new_d$V1), max(new_d$V1), length.out = 50),
 PC1 = mean(new_d$PC1),
 PC2 = mean(new_d$PC2),
 PC3 = mean(new_d$PC3))
new_data$predicted <- predict(lmod, newdata = new_data) # предсказанные значения
new_data$SE <- predict(lmod, newdata = new_data, se.fit = TRUE)$se.fit # стандартные ошибки
new_data$upper <- new_data$predicted + 2 * new_data$SE # верхняя граница доверительной области
new_data$lower <- new_data$predicted - 2 * new_data$SE # нижняя граница доверительной области

Визуализируем влияние первой компоненты.

ggplot(data = new_data, aes(x = ERBB4_N, y = predicted)) +
 geom_ribbon(data = new_data, aes(ymin = lower, ymax = upper), alpha = 0.3) +
  geom_line(color = 'blue', size = 1) +
 geom_point(data = new_d, aes(x = ERBB4_N, y = PC1)) +
 labs(y = 'PC1', x = 'ERBB4_N')

Построить модель не получилось, облако точек явно меняется по другой закономерности.

Задание 4. Построить РСА.

m <- length(colnames(d))

# уберем колонки c не белками
prot_only <- as.matrix(d[,c(-1, -(m-3):-m)])
# стандартизируем переменные
# Function to standardize numeric variables
z_numeric_vars<-function(data_frame) {
  for (n in names(data_frame)) { 
    if (class(data_frame[[n]]) == "numeric" | class(data_frame[[n]]) == "integer"){
      var = paste(n,"_z", sep="")
      data_frame[[var]] <- scale(data_frame[[n]], center = TRUE, scale = TRUE)
      data_frame[[n]] = NULL
    }
  }
  data_frame
}
prot_st <- z_numeric_vars(prot_only)
# заменим NA на среднее
prot <- na.aggregate(prot_st)